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Abstract 

We present a multiscale integrator for Hamiltonian systems with slowly varying quadratic 
stiff potentials that uses coarse timesteps (analogous to what the impulse method uses 
for constant quadratic stiff potentials). This method is based on the highly- non-trivial 
introduction of two efficient symplectic schemes for exponentiations of matrices that only 
require 0(n) matrix multiplications operations at each coarse time step for a preset small 
number n. The proposed integrator is shown to be (i) uniformly convergent on positions; (ii) 
symplectic in both slow and fast variables; (iii) well adapted to high dimensional systems. 
Our framework also provides a general method for iteratively exponentiating a slowly varying 
sequence of (possibly high dimensional) matrices in an efficient way. 



1 Introduction 

One objective of this paper is to obtain an explicit and efficient numerical integration algorithm 
for the following multiscale Hamiltonian system: 
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(1) 



where q slow 



slow 



and qf ast ^pf ast are slow and fast degrees of freedom (in the sense that slow 
degrees of freedom have bounded time derivatives, whereas time derivatives of fast ones may grow 
unboundedly as e — > 0). Observe that a direct numerical integration of ([lj becomes prohibitive as 
e I 0. Notice also that not all stiff Hamiltonian systems are multiscale, and whether a separation 
of timescales exists depends on specific forms of V(-), U (•) and initial conditions. To the authors' 
knowledge, a generic theory that determines whether a stiff system is multiscale has not been 
fully developed yet. 

We will mainly discuss and analyze the case where U (q^ :<lst , q slow ) = ^[qf ast ] T K (q slow )qf ast , 
which we call a quasi-quadratic potential throughout this paper. In this case, the proposed 
method will be able to integrate the system using a coarse timestep. Notice that if K remains 
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constant with respect to q slow , then the impulse method [T5J [3^1 HH 133] allows for an accurate 
and symplectic (see for instance [TB] for a definition) integration of (Q]) using coarse steps. The 
impulse method can, in principle, integrate the situation where if is a regular function of slow 
variables; however, its practical implementation requires a numerical approximation to the stiff 
system 



which generally needs to be based on a numerical integration with small steps. The advantage 
of the impulse method over Verlet is that VV only needs to be evaluated at coarse timesteps, 
but nevertheless its computational cost blows up as e —> 0. 

To use a coarse integration timestep independent of e, we adopt a splitting approach to treat 
the slow and fast variables separately. At each coarse step, we will require an exact solution or 
a numerical approximation to the following stiff system: 



in which q slow is fixed (different from the impulse method). To obtain such an approximation, 
we compute the exponential of a matrix that depends on K(q slow ) and dK(q slow ). Still, if 
not handled appropriately, the cost of this method blows up rapidly as e decreases and/or the 
dimension of the system increases. Furthermore, symplecticity would also be jeopardized by 
inaccuracies of the numerical exponentiations. 

In this paper, we propose an integrator well-adapted to high-dimensional systems, which 
computes the exponentiation in an efficient and symplectic way. Only 0(n) matrix multiplication 
operations at each coarse time step are needed, where n is a preset small integer at most loge" 1 . 
Although simple in appearance, to guarantee the symplecticity (in all variables) of the resulting 
method without diagonalizing K(q slow ) (which is expensive) is a surprisingly difficult problem, 
and in fact it is highly non-trivial even when K(q slow ) is a scalar [22] . 

In addition to a solution to this problem, this paper also provides a general method for 
iteratively exponentiating a slowly varying sequence of (possibly high dimensional) matrices in 
an efficient way (see Section 12.41 and the Appendix) . This method works for any matrices, and 
it is not restricted to the integration of ([I). The preservation of symplecticity associated with 
these two proposed matrix exponentiation schemes (both suit high dimensional systems; the first 
one is in Section [23]) is a core difficulty addressed in this paper. 

Also, useful discrete geometric structures regarding the flow map of a parameter-dependent 
vectorial harmonic oscillator have been studied; for instance, derivatives of its flow map with 
respect to the parameter can be computed using matrix exponentials. 

Although backward error analysis (relating symplecticity and energy conservation) does not 
apply directly to stiff systems (due to large Lipschitz constants), we numerically observe im- 
proved long time behaviors for the proposed integrator, such as near-preservation of energy and 
conservation of momentum maps. We also note that modulated Fourier expansion [5] has been 
proposed to explain favorable long time energy behaviors of some integrators for oscillatory 
Hamiltonian systems. 

We prove the uniform (in e" 1 ) convergence of the method and bound the global error on 
qslow ky (jjj ^ where H is a coarse integration timestep and C is a constant independent of er 1 . 
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2 The proposed method 



2.1 General methodology 

Consider a Hamiltonian system with the following Hamiltonian: 

H(q,p) = \ V J 'M-'p + Viq) + e- x U{qf as \q slow ), (4) 

where < e < 1, q G R d , p G R d . q fast G R df and q slow G R d " are fast and slow variables. 
For the sake of clarity, we will assume that q = (q' ast , q slow ) (and hence d = df + d s ), but the 
method presented here can be generalized to the situation where (q-> ast , q slow ) = 77(g), where ij 
is a diffcomorphism explicitly known beforehand (for instance, see [7] for an example where a 
diffcomorphism puts an extensible pendulum into a quasi-quadratic form) . The idea of separating 
the variables is that fast variables need to be integrated using o(y / e) timesteps by a single scale 
integrator, whereas slow variables can be resolved using o(l) steps. Dividing variables into 
different timescales is a widely used technique in studying stiff and multiscale systems (e.g., 
[8j[T9]). A rigorous definition of separation of timescales can be found in [32], for instance. 

Without loss of generality, we can further assume M to be the identity matrix. The governing 
ODE system is 



'fast ~,fast 

q J — p J 

■fast _ ,-1 dU dV 

P — t dqj ast dq fa S , 

qSlow pSlow 

pslow _ e -l dU dV_ 



(5) 
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which can be split into a sum of three vector fields: 

( qfast _ pfast 
■fast _ _ c -l_dU 
P - fc Q q la.st 

qslow = Q 
■slow _ ,-1 dU 
p — t dqslou , 

such that the exact flow of each could be obtained, which is also symplcctic at the same time. 
Indeed, denote the flow maps of all systems by <p' (s),i = 1.2,3 over a time of s. It is easy to see 
that they are all symplectic. 

Observe that (j) 1 and </> 2 are analytically available. We only consider the case where 3 is 
also analytically or numerically known; more precisely, the numerical solution </r has to have a 
consistent uniform local error over a coarse time step H — o(l), i.e., \\(f) 3 (H) — 3 (i/)|| < CH 2 
for a constant C independent of e^ 1 . This can be satisfied for arbitrary U(-) by a symplectic 
integration with a microscopic timestep h = o(y / e), which is in the same spirit as the impulse 
method. One step update of the proposed method is obtained by composing the three flow maps: 
4> 1 {H) o 4> 2 {H) o <fi 3 (H). Notice that any split can result in a convergent numerical scheme, but 
this particular split treats two timescales independently and therefore is uniformly convergent at 
least in the quasi-quadratic stiff potential case (illustrated later); also, it results in a symplectic 
scheme. 



Remark 2.1. If there were no slow variable, we would compose the flows of ■ 



qfast _ pfast 
•/asi r — 1 dU 



\q fast =0 . 
and < . j ast Q V and obtain a first-order version of the original impulse method. 



3 



Remark 2.2. There are also alternative higher-order ways of composing these flow maps; see, 
for instance, J 161 I #6) /. In fact, the original impulse method is second-order and can be con- 
structed from a second-order composition scheme. However, we will stick to first-order Lie- 
Trotter (4> 1 {H) o <fi 2 (H) o (f) 5 (H) ) in this paper. 



2.2 Quasi-quadratic fast potentials 

We will, from now on, discuss and analyze an analytically solvable case (so that a uniform coarse 
timestep could be used), in which U — h[qf ast ] K{q slow )qf ast , where K is a positive definite 
df-by-df matrix valued function. This fast potential represents stiff harmonic oscillators with 
non-constant but slowly varying frequencies. We call such potentials quasi-quadratic. In this 
case, the first split vector field is 



•fast 

1 

xfast 



fast 



V 

_ e -lf(fqSlow\qfa,st 



= 



(6) 



= — t~ l j[qf ast ] v K{q slow )qf ast 



where the last equation is understood as p. 



slow 



-1 1 



[qf ast ] d l K{q slow )qf ast for i = 1, 



, d,. 



The flow of this dynamical system on qf ast and pf ast is just an exponential map, which in 
this case corresponds to linear combinations of initial conditions with trigonometric coefficients. 
For p slow , because q slow (and hence VK{q slow )) is fixed, one could obtain its exact flow by 
analytically integrating a quadratic function of trigonometric functions. 

When df = 1, the exact flow map of ^ over time H is (letting u> = y e^ 1 K(q slow )): 



jfast 
n fast 



~slow 



i — ^ cos(u;iJ)(^ ast + sin(wTJ) /cop fast 
i y —ujs\x\{{jjH)qf ast + cos{uH)pf ast 
i — ^ q 



slow 



l VK(q slow ) J ^(2u;(H[pf a 

%st n fast 



+ pf ast qf ast 4- to 2 H[qf ast ] 2 ) 



-2up fast q fast cos(2loH) + (-[p 



_\ n fast]2 



+ u z [q 



2\ n fast]2 



) sm(2uH)) 



(7) 



where again the last equation is understood as 



„slow 



H> p . 



slow 



st]2 



pfastqfast 



uj 2 H[qf ast } 2 ) 



i ^d l K(q slow )j^(2oj(H[pf a ' 
-2u}p fast q fast cos(2ujH) + {-\pf ast f + u 2 [qf ast ] 2 ) sm{2uH)) 



(8) 



When df > 2, the obvious method to obtain the exact flow of © is based on a diagonalization 
ofK. More precisely, since if is symmetric, we can write e- 1 K{q slow ) = e- 1 Q(q slow ) T D(q slow )Q(q slow ), 
where e- 1 D{q slow ) = diag[w?, . . ]). Then 



exp 



HI 




\Q T o " 


exp { 


e^HKiq 31 ™) 


)" 


Q T 





HI' 




Q 0" 




\Q T 


" 


e~ 1 HD 


) 


Q 










diag[cos(wiii), . 
diag[— sin(wiiJ)wi, . 



. ,cos(uj df H)} diag[sin(cJiff)/cJi,. 
,-sm(u}d f H)wdf] diag[cos(wii?), 



. ,sm(u} df H)/uj dj ] 
. ,cos(uj df H)) 



(9) 

A similar (but lengthy) calculation will give the expression of the flow on p slow . 

If the diagonalization frame of K(-) is constant, i.e., Q does not depend on q slow ; then Q 
needs to be computed only once throughout the simulation, and then the calculation of the flow 
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on qf ast and pf ast is dominated by the cost of 2 matrix multiplication operations per coarse step 
(at expense per multiplication by the state-of-art Coppersmith- Winograd algorithm 

[5]). However, if the frame varies (Q depends on q slow ), then diagonalizing K at each time step 
can offset the gain obtained by the macro-time-stepping of the algorithm. This is especially true 
if df is large. Moreover, errors in numerical diagonalizations may accumulate and deteriorate 
the symplccticity of cf> 3 . 

In this paper, we address those difficulties by proposing a method, described below, for the 
numerical integration of ((SJ that is symplectic and that remains computationally tractable in 
high-dimensional cases (large df). 



2.3 Fast numerical matrix exponentiation for the symplectic integra- 
tion of © 

The proposed method is based on matrix exponentiation. We will first describe its analytical 
formulation, and then present an accurate numerical approximation that is both symplectic and 
computationally cheap. 

The first step of our method is based on the following property of matrix exponentials illus- 
trated in [2 7) : if N and M are constant square matrices of the same dimension, then 



exp 



-N T M 
N 



H = 



F 2 (H) G 2 {H) 
F 3 (H) 



ith 



F 2 (H) = cxp(-N T H) 

F 3 (H) = exp(NH) 

F 3 {H) T G 2 {H) = exp(N T s)Mexp(Ns)ds 

Therefore, ordering coordinates as qf ast ^pf ast ^ taking N := 



(10) 



(11) 



e- 1 d i K(q slow ) 

obtain that if 



and Mi 



-e- x K{q slow ) 

with i = 1 , . . . , d s which indicates the component of the slow variable, we 



F 2 {H) G 2 ,i(H) 
F 3 (H) 



exp 



N 



H 



then the (linear) flow map on q' ast , pf ast is given by 

exp(NH) = F 3 {H) 

and the drift on p slow is given by 

rt+H i-H r f„.f/.%n T 



1 fast {s) T e' 1 d % K{q slow )q fast {s) ds = 
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p fast (t)_ 
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(12) 
(13) 

ds 
(14) 
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Therefore, 4> 3 (H) is given by: 
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where in the last equation i = 1, . . . , d s . 

In addition, our specific choice of Mi is a symmetric matrix for each i, because K{-) is 
symmetric. Consequently, exp(N T s)Mi exp(iVs) is symmetric, and therefore 



F 3 (H) T G 2 4H) = {F 3 (H) T G 2 ^H)) T 



(16) 



Assuming we have F3 and G 2 ,i (which will be given by Integrator^, Q can be integrated 
by the following: 

Integrator 1. Symplectic multi-scale integrator for (|4|) with U = ^[q^ ast ] T K(q slow )q :i:ast . Its 
one-step update mapping qk,Pk onto qk+i,Pk+i with a coarse timestep H is given by: 
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where F 2 ,k, G 2i k.i (i = 1, • • • , d s ) and F 3t k are numerical approximations of that in (1121) at each 
time step k' (using q k d , ow ), for instance computed by Integrator^ 

To numerically approximate the above flow map (1151) . i.e., to obtain F 3t k and G 2t k,i, we need 
to ensure two points: (i) an approximation of the matrix exponential (and hence F 3 ^k and G 2i k,i) 
will not affect the symplecticity of the resulting approximation of 3 ; (ii) the numerical compu- 
tation of the exponential (fT"2")) will not off-set the savings gained by using a coarse timestep. It 
is highly non-trivial to satisfy both simultaneously, because most matrix exponentiation meth- 
ods will ruin the symplecticity of c/> 3 unless high precision (much higher than the requirement 
on accuracy) is enforced, but then the computational cost will be high. In fact, a necessary 
and sufficient condition for symplecticity is given by Lemma 13.11 and it is unclear how most 
matrix exponentiation methods, for instance those based on matrix decompositions (e.g., diago- 
nalization, QR decomposition) with computational costs of Cd 3 flops, will satisfy this condition 
(unless C is very large and the approximation is very accurate). Also, here is an illustration of a 
popular non-decomposition-based exponentiation method that fails to satisfy this symplecticity 
condition: 



G 



Example: MATLAB function 'expm' [T7] uses a scaling and squaring strategy based on the 
following identity: 

exp(A) = [exp(A/2™)] 2 " (19) 

where n is a big enough preset integer such that X/2 n has a small norm, and therefore Pade 
approximation [17] could be employed to approximate exp(A/2"). The simplest (1,0) Pade 
approximation, which is essentially Taylor expansion to lst-order, gives 



exp(A) a [I + X/2 71 ]' 



(20) 



However, this approximation is not symplectic. For instance, consider a counterexample of 
Obviously, this corresponds to a vectorial harmonic oscillator, and exp(A) 



X 



fl 2 

ought to be symplectic. However, it can be easily checked that A := I + A/2" does not satisfy 
A T J A = J and hence is not symplectic. □ 

Our idea is to obtain F%j. and F^ k using a modified scaling and squaring strategy, in which 
the Pade approximation is replaced by a symplectic approximation originated from a reversible 
symplectic integrator (we use Velocity- Verlet). More precisely, suppose h > is a small constant, 
then we have the following identity: 

H/h 



F2,k(H) G2,k,i(H) 
F 3tk (H) 

F3 k(h>) can be approximated by the following: 



F 2 ,k{h) G 2 ,k,i(h) 
F 3 , k (h) . 
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(21) 



(22) 



which can be easily checked to be symplectic thanks to the specific 0(h 2 ) and 0(h 3 ) corrections 
in the above expression. 

It is a classical result (global error bound of Velocity- Verlet) that links F 3t k{H) with the 
approximated F 3 ^{h): 
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for some constant C > 0, because the approximation in 
Velocity- Verlet integrator with updating rule: 



< e" 1 C*exp(C7?)/i 2 
(23) 



< 



- x l+ i + 



Xi 



corresponds to the celebrated 



(24) 



for the system 



which is well-known to have a 2nd-order global error. 



x = V 

y = -e- 1 K{q s k l , ow )x 

We can repeat the same procedure to get an approximation of F2,k{H) by using the following 
approximated F2 } k(h)' 



exp 





-hi 



he^^iqf, ^) 
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(25) 
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To approximate G2,k,i{h), we follow the result of Lemma 13.11 that in the continuous case 
G%k,i = - J 3,i. f M and let 



G 2 ,k,i{h) = -Jd t F 3tk (h) 



h 1 3 rs(„slow\ 



(26) 



4 - 

Notice that if (1,0) Pade approximation (i.e., lst-order Taylor expansion) is used, we will get 











(27) 



Naturally, (121)1) is a higher order correction of this. 

G2.k.i{H) will also be accurate: since the accuracy of (l20l is well established, the higher order 
corrections that we add in F2k(H), F^ k {H"), G2. k .i(H) will not lead a scheme less accurate. This 
can immediately be seen in the context of the numerical integration of a stable system, where 
a local error of 0(h 2 ) will only lead to a global error of at most e CHh [3T]. We also refer to 
Appendix A in [53] for an analogous error analysis if one prefers to directly work with matrices. 

To sum up, the following numerical approximation of F% yk and G 2 ,k,i will simultaneously 
guarantee symplecticity, accuracy, and efficiency: 

Integrator 2. Matrix exponentiation scheme that complements the updating rule of Integrator 
[7J n > 1 is an integer controlling the accuracy of the approximation of the matrix exponentials, 
k is the same index as the one used in Integrator^ and the following needs to be done for each 
k: 



1. Evaluate K k := K(q s k l , ow ) and diK k 

A k := 



and for i = 1 



^^K{q s k l r ). Let h = H/T 



, . . . , d s , 



hil-e- 1 ^ 



e- x K k h 



Cfc :— 



-e- X K k h I-er l K k \ 



(28) 
(29) 



(30) 
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C k , then repetitively apply 
for j = 1, . . . ,n. 



LT 2,fc,i- r 3,fc 
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1 - T 2,fc,i 
pj + l 
r 3,k 



3. Define F 2>k := F?+\ G 2>k , 



/-ro+1 p c 1 



n+1 
3,fe • 
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Remark 2.3. The trick for the computational save is that raising to the 2 n th power is computed 
by n self multiplications, which is due to the semi-group property of the exponentiation operation. 
An obvious upper bound to guarantee accuracy is n < Cloge -1 (because the error of numerical 
exponentiation is bounded by e~ 1 Ch = e~ 1 CH/2 n ). In all numerical experiments in this paper, 
n = 10 worked well, which is a value much smaller than loge -1 , and this choice of n makes 
the computation cost of the same order as if K could be diagonalized by a constant orthogonal 
matrix. 

Remark 2.4. Observe that, for a finite-time simulation, the cost of computing 4> 3 numerically 
with microscopic time-steps blows up with a speed o/0(e _1 ), whereas the cost of matrix expo- 
nentiations via Integrator^ blows up at a maximum speed o/C(loge _1 ). 

Theorem 13.11 shows that Integrator [5] not only ensures i^.jfc and Fs t k to be symplectic, but 
also guarantees a symplectic approximation to <ft 3 (Eq. I15[) . 

Speed-up is obtained because at each step the computation cost is dominated by 2(d s + 
l)n matrix production operations (of df x df matrices), where n is a small integer. If the 
Coppersmith- Winograd algorithm is used to realize the matrix multiplication operation, then 
the time complexity for exponentiation at each step is nO(d 2 f 376 ) (assuming d s = 0(1); the 
problem of matrix exponentiation is less difficult otherwise) . 

2.4 Fast numerical matrix exponentiation for the symplectic integra- 
tion of ()6]): an alternative 

An alternative way to approximate the flow map (|15p is to use the slowly varying property of 
K to generate a symplectic update of the exponential computed at the previous step. The main 
idea of the method is as follows: given a sequence of matrices {X k } that vary slowly, use the 
approximation 

exp(A fc ) = [exp(AV2")] 2 " « [exp(X fc _ 1 /2") exp((A fc - AVi)/2")] 2 " (31) 

where n is a preset constant. Again, we use the trick of self- multiplication for computing the 
2 n th power, and efficiency is guaranteed exactly as before. 

Accuracy is achieved because, as shown in the following theorem, the approximation error 
decreases at an exponential rate with respect to n. 

Theorem 2.1. Theorem 5 in '251: 

\\exp(A + B) - (exp(A/2 Tl )exp(B/2™)) 2 " || 2 < 2-"- 1 e max ^ A+s )<^ A > + ^ s » || [A, B] || 2 (32) 

where fJ,(X) is the maximum eigenvalue of (X* + X)/2, and [A, B] = AB — BA is the canonical 
Lie bracket. 

Remark 2.5 (Generality). This exponentiation method based on corrections pip is not limited 
to the integration of but works for repetitive exponentiations of any slowly varying matrix. It 
would also work for a set of matrices, as long as they could be indexed to ensure a slow variation. 

For our case, X k and A are identified with N in Section 12.31 at each timestep, and B is 
identified as the difference in iV's between consecutive steps. Since K(q slow ) (and hence N as 
well) is changing slowly, ||B||2 *C ||^4||2i furthermore, the calculation of [A,B] (omitted; notice 
that B is nilpotent) shows that ||[A,-B]|| -C Therefore, the error bound here (1321) is much 

smaller than that based on scaling and squaring for the same n. Consequently, we will be able to 
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further decrease the value of n by a few (not a lot because a decrease in n exponentially increases 
the error). 

-N T Mi 



The reason that we do not identify Xk and A with 







N 



is due to a consideration of 



symplecticity in all variables, because otherwise G% t k,ii obtained as the upper- right block of the 
exponential, will not be exactly the derivative of i*3,fc. Instead, we let C?2 fe j = ~ J a ircs F- 



3,k, 



where F^k is updated from F^^-i using (pH) . Taking the derivative, however, incurs addi- 
tional computation, because now depends on not only qf. but also qf. 1 "™, and therefore 
®Qk°r / ®Qfk-i)' j nas to ^ e com P u ted so that a chain rule applies to facilitate the computation. 



In the end, the computational saving based on updating the exponentiation becomes less signif- 
icant due to the extra cost in updating dqf},°™/dq 
convoluted. We leave the details to the appendix 



icant due to the extra cost in updating dq^f /dqf^U, p but the implementation becomes more 



3 Analysis 

3.1 Symplecticity 

For a concise writing, we carry out matrix analysis in block forms in this section. Coordinates are 

J 0" 



ordered as qf ast 7 pf ast , q slow ^p slow ^ anc } therefore 



J 



is the coordinate representation of 



the canonical symplectic 2-form on the full phase space (abusing notations, we use J := 



/ 
-I 



to represent the symplectic 2-form on both the fast subspace (for qf ast 7 pf ast ) and the slow 
subspace (for q slow 7 p slow )- 7 this should not affect the clarity of the analysis). We also recall that 
a map x n- <fi(x) is symplectic if and only if (f>'(x) T $cf)'(x) = JJ or cf>' (x) T J(j>' (x) = J for all x's 
(depending on whether x represents all variables or only slow or fast variables). 



Lemma 3.1. The numerical approximation to <p given by (|18p is symplectic on all variables if 
and only if F$ t k is symplectic and, for i = 1, . . . , d s , G2,fe,i = — J Jpril; (note that for a fixed i, 



G 



2,kA, 



and J are df x df matrices). 



Proof. For conciseness and convenient reading, write q{.? st and p{.? st as qf and pf, d/dqf. l ,°™ as 
di, and G2,k,i and F^^ as G2A and F3 in this proof. 

The Jacobian of the numerical approximation to 
be computed as: 



fast 



qk',Pk' i-> qk+i,Pk+i given by (O can 



.4 




diF 3 



d d F, 



I 





— * 





(33) 



^0 0y 
{qf ,> } )F!G 2 , 

where = ^[qf,pf] T dj(F^G2,i)[qf,Pf], and the 0's in the upper right block, the lower left 

block, and the lower right block respectively corresponds to d/-by-l, 1-by-d/, and d s -hy-d s zero 
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matrices. Notice that we have F^G 2 ,i in the lower left block because F^G 2 ,i is symmetric (their 
exact values satisfy this because of (|16[) . and their numerical approximations satisfy this because 
of Lemma 13.61) . 

Symplecticity is equivalent to A T ]sA = J, whose left hand side writes out to be 



A JA 



F?J 




Pf) diF% J 


iff 


P T f)d d FjJ 




fO 0\ 




{o o) 



G 2,1^3 



G 2,d s F 3 



T 

* 


I 


-I 






x A 



F[JF 3 + 


{FiJd 1 F, + Gl 1 F,)^ ... {FfJd d 


S F 3 + G% ds 


F 3 ) 


GO 


(o 




A 


{'If-Pf 1 <).!..! Ji),l :i <!,: pr ). . 


d s ;j=l,...,d 3 





— * 
H 


T +* 


i 










-J 






(34) 

where A is naturally negative the transpose of the upper-right block because A T SA is skew- 
symmetric for any A. 

This is equal to J if and only if the upper-left block and the bottom-right block arc both J and 
the upper-right block and the bottom-left block are both zero. The requirement on upper-left 
block is 

Fj JF 3 = J (35) 

By the arbitrariness of qf and pf, the requirement on upper-right and bottom-left blocks trans- 
lates to: 

F 3 T J^F 3 + G^Fs = (36) 

which further simplifies to 

G 2 ,i = -JdiF 3 (37) 

because F 3 T Jd l F 3 = d l (F 3 T JF 3 ) - d t F^JF 3 = -d l F^ T JF 3 , F 3 is invertible due to {35]), and 
J T = -J. 

The bottom-right block needs to be J, and this requirement is equivalent to 



, ir .p r \ ' ( diFfjdjFs + l -d l {F^G 2 . J ) - ^(if G 2ti ) ) [q f ;p f ] = 



(38) 



By (I3T1) . the above left hand side rewrites as 
[q f -p f ] T (diFTjdjFs - ^F^JdjFa - ^Jd^F, + ^if Jc^ + ^Jd^F^j [q } ;p f ] 

= tef,Pf} T (^ifj^-Fa) [q f] p f ] + [q f ;p f ] T Q^ifJc^) [q J]Pf ] (39) 

Since what are summed up above are just two real numbers, the second number remains the 
same after taking its transpose, which due to J T = — J yields 

ilflPff (\d 3 FiJd t F 3 \ [ qf ; Pf ] = -[qf,p f ] T (\d^ ' Jd,F 3 \ [q f ;p f ] (40) 
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Therefore, 



does hold. □ 



□ 



Lemma 3.2. In Integrator^ all A k and C'k are symplectic; moreover, all F^^k and F 3 ^ are 
symplectic, too. 

Proof. Straightforward computation using (|2"5)) and (|2T))) shows that A\J A k = J and C k JCk = 
J. Moreover, since the product of symplectic matrices is symplectic, all i^.fc and F 3> k, being 
powers of Ak and Cfc, are symplectic. □ □ 

Lemma 3.3. In Integrator^ A\Ck = I (and equivalently C k Aj = I) for all k: moreover, 
F^kFztk = I (and equivalently Fz k^ik = 

Proof. Straightforward computation using (|28p and Q29]) shows that AjCk = I. Therefore, 
{AkAkfCkCk = A T k IC k = I, and by induction (A 2 k ") T Cf = I, i.e., F? k F 3 , k = I. □ 



Lemma 3.4. In Integrator^ Bk,% = - 
all k and i. 

Proof. Use the short-hand notation d% 
O shows that B kfi -- 

F2.k G2,k, 

F 3 , k 



d 



Since 



JdiCk for all k and i. 

A k Ik-. 1 

o c k 



C k for all k and i, and G 2 ,k,i = ~J d sio-^ F^k for 

k 1 ,i 

— . Straightforward computation using (|3"0|) and 



for all i, by induction, it is only necessary to prove that 
A k B k ;+B k iC k and F 3<k = C k C k , and the 



G2,k,i = —JdiF^^k when n = 1. In this case, Gi.k 
equality can be proved by the following: 

Because Bk.i — —JdiCk, C k Ak = / (Lemma 13. 3[) and J = C k JCk (Lemma 13. 2[) . we have 

G T k A k B Kl = -ClJCkd % C k (41) 

Since symplectic matrix is nonsingular, this is 

AkBk : i = —JCkdiCk (42) 

Adding Bk^C k = -Jd l C k C k , we have 

AkBk,i + Bk,iCk — —Jdi(CkCk) (43) 

Hence, the induction works. □ □ 

Lemma 3.5. In Integrator^ C k Bk.i — B ki Ck for all k and i. 

Proof. This can be shown by straightforward computation using (|30l) and (1291) . □ 
Lemma 3.6. In Integrator^ F 3k G2, k ,i — G\ k iF 3 ,k for all k and i. 

Proof. By Lemma 13.51 C k Bk,i = B^-Ck for all k and i. By Lemma 13.31 AjCk — I and 



CTA k 



I. 



Since 



F2.k G2,k,i 




'A k 


Bk,i 


. F 3,fe _ 










2" 



for all i, by induction, it is only necessary to prove that 

F 3 k G2,k,i = G\ k i F z .k when n = 1. In this case, G 2 , k ,i = A k B kii + B k<i C k and F 3 . k = C k C kl 
and this equality can be proved upon observing for all i: 

GkGk(AkBk t i + Bk,iCk) — C k B k ,i + C k C k B k ,iC k = B ki C k + C k B k i C k Ck 

= Bk,iAkCkCk + B k { CkCk = {AkBkj + Bk,iCk) T CkCk (44) 



□ 



12 



Theorem 3.1. The proposed method (Integrators^) is symplectic on all variables. 

Proof. By Lemma 13.21 13.41 13.11 and 13.61 the numerical approximation to cf> 3 given by (fT5)l is 
symplectic on all variables. 

The flow given by (fTT)) is symplectic on all variables as well, because it is the composition of (p 1 
and (j) 2 , which respectively correspond to Hamiltonians r H\{q^ ast ,p fast , q slow ^p slow ^ — [p stow ] 2 /2 
and U2(q fast ,p fast ,q slow ,p slow ) = V(q fast , q slow ), and hence both are symplectic. 

Consequently, the proposed method, which composes (IT7|) and (fT5|) . is symplectic. ■ □ 



3.2 Uniform convergence 

This integrator is convergent due to splitting theory [35], i.e., the global error on q slow ; qi :<lst ) p slow ] pi »«* 
is bounded by e~ 1 CH for some constant C > in Euclidean norm. 

Moreover, this integrator is uniformly convergent in q under typical or reasonable assump- 
tions, and hence H can be chosen independently from e for stable and accurate integration. 

Condition 3.1. We will prove a uniform bound of the global error on position for Integrator]^ 
under the following (classical) conditions: 

1. Regularity: In the integration domain of interest, VV(«) is bounded and Lipschitz continu- 
ous with coefficient L, i.e. ||VV(a) — W(i>)||2 < L\\a — b\\2- 

2. Stability and bounded energy: For a fixed T and t < T, denote by x(t) = (q(t),p(t)) the exact 
solution to l[5]). and by Xt = {qt,Pt) the discrete numerical trajectory given by Integrator^ 
then \\x(t)\\l < C, \\x t \\% < C, \H(q(t),p(t))\ < C and \H(q t ,p t )\ < C for some constant C 



independent of e 1 but dependent on initial condition 



and possibly T as well. 



Condition 3.2 (Slowly varying frequencies:). Consider the solution q(s),p(s) up to time s <= H 
to the system 

dq fast = p fast dt 

dq slow = p slow dt 

dp fast = -dv/dq fast (q fast ,q slow )dt- e~ 1 K(q slow )qf ast dt ' ' 

dp siow = -OV/dq slow (q fast ,q slow )dt - e- 1 \[qf ast ] T \JK{q slow )qf ast dt 

with initial condition q(0),p(0) in the domain of interest that satisfies bounded energy. Assume 
that qf ast can be written as 

df 

Q{t) V e^eadt) cos[^^e~^e l {t) + 0,] (46) 



i=l 

where Q(t) is a slowly varying matrix (i.e., Qij(t) € C 1 ([0,iJ]) and there exists a C independent 
of e _1 such that \\Q(t)\\ < C and \\Q(t)\\ < C for all t G [0,H]), indicating a slowly varying 
diagonalization frame, df is the dimension of the fast variable, are standard vectorial basis 
of M. d f , ai(t)'s are slowly varying amplitudes (in the same sense as for Q(t)), 0i(t)'s are non- 
decreasing and slowly varying in the sense that 0i(t) € C 2 ([0, H]), \&i{t)\ < C, \0i(t)\ < C, and 
C\ < Oi(t) < C'2 for some C > 0, C\ > 0, C2 > independent of e _1 , and ifn's are such that 
0i(O) = 0. 
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Remark 3.1. In the case of constant frequencies (K(-) being a constant) and no slow drift 
(V(-) being a constant), we have qf ast = QX)f=i ^iV^ a i cos[\/e™ 1 Wii + 4>i] (the amplitude is 
0(yfe) because of bounded energy). When K is not a constant, Condition \ 3.2\ is supported by 
an asymptotic expansion of qf ast . In particular, to the leading order in e, we have 9i(t) — u>i(t) 
where the uf(t) are the eigenvalues of K(qf ow ). The rigorous justification of this asymptotic 
expansion for df > 1 is beyond the scope of this paper. 

Lemma 3.7. // Condition 1 3. £1 holds, there exists C\ > 0,6*2 > independent of e~ 1 such that 
" H / \ 

(47) 



f(t)q fast (t)dt\\ < e Ci max ||/(s)|| + C 2 H max \\f(s)\\+0(H 2 

n - y 0<s<H a<s<H 

for arbitrary matrix valued function f £ C^QO, H]) that satisfies /(0) = 0. 

Proof. Recall the form of qf ast in Condition 13.21 It is sufficient to prove that for all i's the i-th 
component of q' ast satisfies (14T1) , whereas the i-th. component writes as: 

d f 

V~^Qij{t)aj{t) cos[V^6i{t) + (Pi] (48) 

3=1 

Furthermore, since summation commutes with integral and therefore will only introduce a factor 
of df on the bound, it is sufficient to prove (|4"7) for qf ast = yfeQij(t)aj(t) cos[V e~ 1 9i(t) + <f>j\. 
On this token, we could assume that we are in the ID case and absorb Q{t) into cij(t). 

Similarly, slowly varying ai(t) can be absorbed into the test function f(s), and doing so will 
only change the constants on the right hand side. Therefore, it will be sufficient to prove that: 



H 



\/icos[V / F T 6i(t) + <p\f{i)dt 



<e Ci max \f(s)\+C 2 H max \f'(s)\+0(H 2 

1 0<s<H 0<s<H 



(49) 



for a scalar valued function / € C 1 ([0,_ff]) that satisfies /(0) = 0. 

By Condition 13. 2[ 9 is strictly increasing. If we write r = 9(t), there will be a 9~ 1 such that 
t = 9~ 1 (t). With time transformed to the new variable r, the integral on the left hand side of 
(gSJl is equal to 

V~ecos[V^T + 0/(0-i( T ))2l_(r) dr 



By integration by parts, this is (since /(0) = 0) 



esinlVe^H + <f>]f(H) 



9(H) 



8(5) 



df (d9~ x 
dt \ dr 



f{0-\r)) 



d 2 9- 
dr 2 



Because 9 < C ', lo - CH < 9 < w + CH, where w := (9(0) > C x > 0. Together with ^ 
have ^t— =l/u + 0(H). Similarly, we also have 

d 2 9- 1 d 1 dt d 1 
dr 2 ~ dr 0(f) ~ dr dt 0(f) 



= ^-f -A- = —J—9(t) = Oil) 
9(tf 



(50) 



-(r) 



(51) 



(52) 



It is easy to show that 9(H) = 0(H). Together with sin(-) being 0(1), the left hand side in ([4"9"]l 
is bounded by 



ef(H)0(l)+eO(H)[0(l) max J/(s)| + 0(1) max \f(s)\ < e 0(1) max \f(s)\ + 0(H) max |/(s)| 



0<s<ff 



0<s<ff 



0<s<H 



0<s<H 

(53) 
□ 
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Theorem 3.2. If Conditions \3.1\ and \3.2\ hold, the proposed method (Integrator^ for system 
(O has a uniform global error of 0{H) in q, given a fixed total simulation time T — NH: 



\q(T)-q T \\ 2 < CH 



(54) 



where q(T),p(T) is the exact solution and qr,PT is the numerical solution; C is a positive 
constant independent of e _1 but dependent on simulation time T, scaleless elasticity matrix K, 

slow potential energy V(-) and initial condition | ^° ||2- 

[Po] 

Proof. Let K be a constant matrix and consider the following system: 
' dqf ast = pf ast dt 

')dt- t- x kq^ ast dt 
•)dt 



dq sl ° 

d ~slow 



_ _Qy I Q^slaw Iq-fast fjslow\ 



(55) 



Integrator [TJ applied to the system ((55)) under Condition 13.11 has been shown in 33J to be 
uniformly convergent in "energy norm" (or equivalcntly, uniformly convergent on position and 
non-uniformly convergent on momentum). Recall that the "energy norm" was defined in [33] to 
be 

\\%p]\\e = 



q T q 



ep T K 1 p, 



(56) 



but in fact K 1 is not important because it is just 0(1), and the following definition would also 
work for the proof there: 

\\[q,p}\\e 



= Vq T q- 



ep T p 



(57) 



Observe that, ([55]) is proportional to the physical energy of ^[i(K~ x l 2 q, K~ 1 l 2 p). 

The system considered here, however, is (1451) . To prove uniform convergence for (|45|) . it is 
sufficient to show that (i) a 5 difference between two trajectories of l|55[) in energy norm leads to 
a difference of 6(1 + CH) in energy norm after a time step H (ii) trajectories of (|5"5"|) and (|4"5")) 
starting at the same point remain at at a distance at most 0(H 2 ) in energy norm after time H, 
i.e., a 2nd order uniform local error, (i) was shown by Lemma 6.5 in (33) . and we will now prove 
(ii). 

We can assume without loss of generality that we start at time 0, and let K = K(q ow (0)), 

qfast,slowfQ\ _ qfast,slow fQ\ ( w h ere gfast,slow _ fgfast ^ gSlow\\ an( j pf ast, slow Iq\ _ pf ast, slow 

We first let x — qt ast — qf ast and y = ftt ast — pf ast j and proceed to bound x and y: 
The evolutions of x and y follow from 



y 



Writing f x = - (^r(q) - gjfi*(qj) and f 2 = (K- K(q slow ))qf ast , we have 




If we let B(t) = cxp 



I 

-e~ x K 



t , we will have 



x(t) 

y(t) 



B(t) 



x(0) 
2/(0) 



B{t 



fi ~ e~ x h 



ds 



(58) 



(59) 



(60) 
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The first term on the right hand side drops off because x(0) = and y(0) — by definition. 

Since K is a constant matrix, it is sufficient to diagonalize it and treat each diagonal element 
individually. Hence, assume without loss of generality that we are in the ID case. Then B(s) = 



i(Ve- l Ks) sin(V 'e~ 1 Ks)/Ve- 1 K 



y(t)- 



As a consequence, 



cosIvV 1 ^ - a)] fi - e~ l {K - K(q slow ))q fast 
By Lipschitz continuity of V (Item [1] of Condition 13. we will have 

\h(t)\<L\x(t)\=L\ f y(s)ds\=0(t) 



ds 



(61) 



(62) 



The first inequality holds because fi is the difference between partial derivatives of V, which 
could be bounded by the difference between full derivatives. The last equality holds because 
y = p — p is bounded due to the fact that [q(s),p(s)] and [q(s),p(s)] are bounded (Item [2] of 
Condition 13. ip . Consequently, we have 



s[Ve^R(t-s)]fids 



< / \fl\=0(t 2 ) 



(63) 



In order to bound J* cos[Vp^K(t-s)} e- l {K - K{q slow ))qf a8t ds. we use Lemma I3T71 (with 

the choice of / = K - K{q slow )). Indeed, cos[V e _1 A"(< - s)] can be absorbed into q fast (s) = 
^/ecos[V e~ 1 9(s) + cj)]: due to an equality 2 cos(^4) cos(i?) = cos(A + B) + cos(A — B), 9 will be just 
added by ± \f~K and <j) will have a new constant value, neither of which will violate Condition 
321 

For /, we clearly have / = at s = 0. By mean value theorem, there is a £ s such that 
f(s) — K o q slow (0) - K o q slow (s) = dKo £'° m (£,) • s, and therefore f(s) = O(s). Similarly, 



f(s) — 0(1). Plotting these two bounds in Lemma \3. 71 we obtain 



cos 



[Vr^(t - s)} le'^K - K{q slow ))q fas 



ds 



<D(t) 



(64) 



Putting this together with (|S3")l . we arrive in y(t) = 0{t), and x(t) — f*y(s)ds — 0(t 2 ) 
follows. 

Next, we bound y: since 



cos 



[\Je- l K{t - s)] \e- l (K - K(q slow ))q fas 



ds 




L 









cos[. . .]e~ 1 0(s)v^O(l) cos[. . .] ds 



(65) 

we have y(t) = e~ 1/2 0(^ 2 ). Together with x(t) = 0{t 2 ), this is equivalent to || [atr, 2/]||js = 0{t 2 ). 



Similarly, we can bound q slow - q slow and p slow - p slow . 



p siow _ pslow ) then we have . 



Let x s = q slow - q slow and y s = 



x = y 
V s =- 



(opL(q) ~ oipL(q)) - e-^ 2 {qf^} T VK(q° l ™)q faS 



(66) 



Analogous to before, the first term on the right hand side of the y s dynamics is 0(t). Since 
qfast _ 0( e i/2^ t ne seconc l term on the right hand side is 0(1). Therefore, y s = 0(1), y(t) = 
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7/(0) + 0(t) = 0(t), and x(t) — x(0) + J * y(s) ds = 0{t 2 ). For our purpose of fast integration, 
we use a big timestep H > y^e, and hence y(H) = O(H) < cT x l 2 0(II 2 ) (notice that if H < \/e, 
we do not even need to prove uniform convergence, because the non-uniform error bound that is 
guaranteed by Lie- Trotter splitting theory is already very small). 

0(H 2 ) and e^ 1 / 2 0{H 2 ) bounds on separations of slow position and slow momentum imply 
a 0(t 2 ) uniform bound in energy norm (analogous to that of the fast degrees of freedom). This 
demonstrates a 2nd-order uniform local error on all variables in energy norm, and therefore 
concludes the proof. □ 

Remark 3.2. Unlike (|54[) . a global bound on the error of momentum will not be uniform. The 
error propagation is quantified in energy norm, and in 2-norm we will only have e _1 / 2 C(i? 2 ) local 
error and e^ x ^ 2 0{H) global error on momentum. In fact, Integrator^ applied to the constant 
frequency system (1551) is non-uniformly convergent on momentum [331. 



4 Numerical Examples 

4.1 The case of a diagonal frequency matrix 

Consider the Hamiltonian example introduced in |22j : 

n = \vU \p\ + (x 2 + y 2 -l) 2 + 1 -(l + x 2 )u 2 y 2 

When uj = eT 1 ! 2 ^> 1, bounded energy translates to initial conditions x(0) « 
satisfy separation of timescales: x is the slow variable, and y is the fast. K(x) = 1 - 



2V1 + 2! 2 



+ 



\/l + X 2 



(67) 

(jjy(0), which 
- x is trivially 

is an adiabatic 



diagonal. In addition to conservation of total energy, / = 
invariant. 

A comparison between Variational Euler and the proposed method is shown in Figure [T] There 
it can be seen that preservations of energy and adiabatic invariant are numerically captured at 
least to a very large timescale. Since there is no overhead spent on matrix exponentiation here, 
an accurate lOOx speed up is achieved by the proposed method (because H/h — 100). 

It is known that the impulse method and its derivatives (such as mollified impulse methods) 
are not stable if the integration step falls in resonance intervals (mollified impulse methods 
have much narrower resonance intervals, which however still exist) j!2l [3]. Similarly, it will be 
very unnatural if the proposed method does not have resonance, because it reduces to a 1st- 
order version of impulse methods when there is no slow variable (Remark I2.1[) . In fact, in our 
numerical investigation (Figure[5]), we clearly observe resonance frequencies before the integration 
step reaches the unstable limit (around H s» 0.5), and widths of resonant intervals increase as 
H grows for this particular example; however, we will not carry out a systematic analysis on 
resonance due to limitation of the length of a short communication. 



4.2 The case of a non-diagonal frequency matrix 

Extend the previous example to a toy example of 3 degrees of freedom: 



1/ l„2 . I 2 . I 2 i / 2 i 2, 2 



2 1X T y 2' 
It is easy to check that eigenvalues of K (x) 



I) 



y 


T 


l + x 2 x 2 - l 




V 


z 




x 2 - l 3x 2 




z 



are both positive when 



l + x z ar - l 
x 2 -l 3x 2 

x > 0.44, which will always be true if the initial condition of x stays close to l and u is big 



(68) 
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The proposed method 



Variational Euler 




0.64 
>> 0.62 
E 1 0.6 
C 0.58 
W 0.56 

0.54 

c 

ra 
to 

> 0.38 

O 0.36 

"fa 0.34 

n3 0.32 











Time 




(a) The proposed method with coarse timestep H 
0.1 



(b) Variational 
0.1/w = 0.001 



Time 



Euler with small timestep h 



The proposed method 




(c) Very long time simulation by the proposed 
method with coarse timestep H = 0.1 



Figure 1: Simulations of a diagonal fast frequency example II67H by the proposed method and Variational Euler. 
w = 100; x(0) = 1.1, 3/(0) = 0.7/u). 



enough. In this case, bounded energy again implies x(0) ~ wy(0) ~ wz(0) and gives well- 
separation of timescalcs: x is the slow variable and y and z are the fast. K (x) has its orthogonal 
frame for diagonalization as well as its eigenvalues slowly varying with time. 

Figure[3]shows a comparison between Variational Euler, the proposed method with the matrix 
exponentiations computed by diagonalization and analytical integration (Eq. [S] diagonalization 
implemented by MATLAB command 'diag'), and the proposed methods based on exponentiations 
(Eq. [10] and [T5)l via MATLAB command 'expm' [17] and via the fast matrix exponentiation 
method (Integrator The default MATLAB matrix multiplication operation is used. All 
implementations of the proposed method are accurate, except that numerical errors in repetitive 
diagonalizations contaminated the symplecticity of the corresponding implementation over a long 
time simulation (as suggested by drifted energy), whereas other two implementations, respectively 
based on accurate but slow 'expm' and fast symplectic exponentiations, do not have this issue. 
In a typical notebook run with MATLAB R2008b, the above four methods respectively spent 
11.12, 0.23, 0.29 and 0.24 seconds on the same integration (till time 50), while 0, 0.14, 0.18, 
and 0.14 seconds were spent on matrix exponentiations. Computational gain by the symplectic 
exponentiation algorithm will be much more significant as the fast dimension becomes higher. 
Notice also that the computational gain by the proposed method over Variational Euler will go 
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Figure 2: Investigation on resonance frequencies of the proposed method on example K71 . The ratio between 
x(T)|t = ioo integrated by the proposed method integration and benchmark provides the ruler: a ratio closer to 
1 means a more accurate integration, and deviations from 1 stand for step lengths that correspond to resonance 
frequencies. Time step H samples from 0.001 to 0.2 with an increment of 0.001. ui = 100; x(0) = 1.1, y(0) = 0.7/oj. 
Benchmark is obtained by fine VE integration with h = 0.01/u;. 



to infinity as e — > 0, even if the fast matrix exponentiation method is not employed. 



4.3 The case of a high-dimensional non-diagonal frequency matrix 

Consider an arbitrarily high-dimensional example: 



u 



+ 7;V T y + {x T x + q 2 - l) 2 + -uj 2 x T T{q)x 



(69) 



where q,p £ M. correspond to the slow variable, x, y £ M. df correspond to fast variables, and T(q) 
is the following Toeplitz matrix valued function: 



T(q) = 



1 


e 


f ■ 


. ff- 1 




i 




. q d f- 2 


f 




1 


q d f-3 


d f -i 


qdf-2 


qdj-3 


1 



(70) 



where q = q/2 so that eigenvectors and eigenvalues vary slowly with q given an initial condition 
of q(0) w 1. Note that the expression of T(-) is highly nonlinear. 

We present in Figure 2] a comparison between Variational Euler and the proposed methods 
with the matrix exponentials computed by MATLAB command 'expm' and by the fast matrix 
exponentiation method (Integrator [2J on a high dimensional example with df = 100. Accuracy- 
wise, the proposed method simulations yield results similar to VE (note that fast variables are not 
fully resolved due to a coarse time step that is larger than their periods). Speed- wise, Variational 
Euler, the proposed methods via 'expm' and via symplectic exponentiation respectively spent 
136.7, 66.0 and 12.0 seconds on the same integration, while 65.7 and 11.7 seconds were spent 
on matrix exponentiation operations in the latter two. Notice that if Coppersmith- Winograd [5] 
is used to replace MATLAB matrix multiplication, the number 11.7 should be further reduced. 
In spite of that, the proposed method with the proposed matrix exponentiation scheme already 
holds a dominant speed advantage, and this advantage will be even more significant if w and/or 
df is further increased (results not shown). 
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Figure 3: Simulations of a non-diagonal fast frequency example (168 6 by Variational Euler, the proposed method 
with different implementations of matrix exponentiations, uj = 100, VE uses h = 0.1/oj = 0.001 and the proposed 
method uses H = 0.1 and n = 10; x(Q) = 1.1, y(0) = 0.2/ui, z(0) = 0.1/cj, and initial momenta are zero. 



5 Related work 

Stiff integration: Many elegant methods have been proposed in the area of stiff Hamiltonian 
integration, and some are closely related to this work. An incomplete list will be discussed here. 

Impulse methods [H2 [35] admit uniform error bounds on positions and can be categorized 
as splitting methods [33]. In their abstract form, impulse methods are not limited to quadratic 
stiff potentials; however, their practical implementation requires an approximation of the flow 
associated with the stiff potential. Our method is based on a generalization of the impulse 
method to (possibly high-dimensional) situations where the stiff potential contains a slowly 
varying component. Although simple in its abstract expression, the practical implementation of 
this generalization (for high-dimensional systems) has required the introduction of a non-trivial 
symplectic matrix exponentiations scheme. 

Impulse methods have been mollified [TU [25] to gain extra stability and accuracy. How- 
ever, mollified impulse methods and other members of the exponential integrator family |14) . 
for instance Gautschi-type integrators |18j . are not based on splitting, and hence the splitting 
approach in this paper does not immediately generalize them. 

The reversible averaging integrator proposed in |23] averages the force on slow variables and 
avoids resonant instabilities. It treats the dynamics of slow and fast variables separately and 
assumes piecewise linear trajectories of the slow variables, both in the same spirit as in our 
proposed method; it is, however, not symplectic, although reversible. 

Implicit methods, for example LIN [57] . work for generic stiff Hamiltonian systems, but 
implicit methods in general fail to capture the effective dynamics of the slow time scale because 
they cannot correctly capture non-Dirac invariant distributions |24) , and they are generally slower 
than explicit methods if comparable step lengths are employed. 

IMEX is a variational integrator for stiff Hamiltonian systems [3T] ■ It works by introducing 
a discrete Lagrangian via trapezoidal approximation of the soft potential and midpoint approx- 
imation of the stiff potential. It is explicit in the case of quadratic fast potential, but is implicit 
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Figure 4: Simulations of a non-diagonal fast frequency high-dimensional example HOI by Variational Euler, 
the proposed method via MATLAB matrix exponentiation 'expm,' and the proposed method via fast matrix 
exponentiations (n = 10). Fast variable dimensionality is df = 100. w = 1000. VE uses h = 0.1/ui and the 
proposed method uses H = 0.1, q(0) = 1.05, x(0) is a df + 1-dimensional vector with independent and identically 
distributed components that are normal random variables with zero mean and variance of 1/oj/ y/df (so that 
energy is bounded), and initial momenta are zero. Only trajectories of the first two fast variables were drawn for 
clarity. 



in the case of quasi-quadratic fast potentials. 

A Hamilton- Jacobi approach is used to derive a homogenization method for multiscale Hamil- 
tonian systems [22,, which works for quasi-quadratic fast potentials with scalar frequency and 
yields a symplectic method. We also refer to [7] for a generalization of this method to systems 
that have either one varying fast frequency or several constant frequencies. The difficulty with 
this elegant analytical approach would be to deal with high-dimensional systems. 

Other generic multiscale methods that integrate the slow dynamics by averaging the effective 
contribution of the fast dynamics include: Heterogeneous Multiscale Methods (HMM) [TTJ [5J SJ 
[U[9], the equation-free method [2Ql H3J HS| , and FLow Averaging integratORS (FLAVORS) [32]. 
Those methods can be applied to a much broader spectrum of problems than considered here. 
However, they all essentially use a mesoscopic timestep, which is usually one or two orders of 
magnitude smaller than the coarse step employed here. Moreover, symplecticity is a big concern. 
In their original form, Heterogeneous Multiscale Methods and equation-free methods are based 
on the averaging of the instantaneous drifts of slow variables, which breaks symplecticity in all 
variables. Reversible and symmetric HMM generalizations have been proposed [2l[29]. FLAVORS 
|32) are based on averaging instantaneous flows by turning on and off stiff coefficients in legacy 
integrators used as black boxes. In particular, they do not require the identification of slow 
variables and inherit the symplecticity and reversibility of the legacy integrators that they are 
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derived from. 

Matrix exponentiation: In the case of quasi-quadratic stiff potentials, the proposed algo- 
rithm exponentiates a slowly varying matrix at each time step. When the elasticity matrix K 
is not diagonalizable by a constant orthogonal transformation, a numerical algebra algorithm is 
employed for that calculation at the expense of 0{n) operations of df-hy-df matrix multiplica- 
tions per time step, where df is the dimension of fast variable (and hence K), and n is a preset 
constant that is at most log(e~ 1 ). 

There are various approaches to exponentiate a matrix, including diagonalization, series meth- 
ods, scaling and squaring, ODE solving, polynomial methods, matrix decomposition methods, 
and splitting, etc., as comprehensively reviewed in |25j . Many of these methods, however, differ 
from our approach here in that they do not guarantee that the resulting implementation of the 
proposed method to be symplectic as it analytically should be, unless high precision (hence slow 
computation) is required; most of them could not even guarantee a symplectic approximation to 
F 2 and F 3 . 

The proposed approach (Integrator [2]) obtains its efficiency by a trick of self- multiplication, 
which was previously used in the method of scaling and squaring j!7j . However, the Pade approx- 
imation used in scaling and squaring is replaced by a symplectic and reversible approximation 
based on the Verlet integrator. Consequently, symplecticity and better efficiency are obtained, 
and accuracy is kept. Improvements by this numerical exponentiation over 'expm' and 'diag,' in 
terms of both accuracy and speed, are observed numerically in Subsection 14.21 and 14.31 

For our alternative approach (see (|3ip for the general strategy and Appendix for implemen- 
tational details for the specific purpose of multiscale integration), it uses the slowly varying 
property of the matrix to repetitively modify the exponential from the previous step by a small 
symplectic change to get a new exponential. Regarding updating matrix exponentials, since 
there are results such as |10) on relationships between perturbed eigenvalues and perturbation 
in the matrix, a natural thought is to use eigenstructures that were explored in the previous 
step as initial conditions in iterative algorithms (such as Jacobi-Davidson for eigenvalues |30] 
or Rayleigh Quotient for extreme eigenvalues [31]). This idea, however, did not significantly 
accelerate the computation as we explored in numerical experiments with an incomplete pool 
of methods. Other matrix decompositions methods (QR for instance) did not gain much from 
previous decompositions either in our numerical investigations. Our way of exponential updating 
is essentially an operator splitting approach, which is analogous to our main vector field splitting 
strategy that yields the proposed multiscale integrator. 
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Appendix: an alternative matrix exponentiation scheme 

We will present in Integrator [3] an alternative (symplectic) way for computing Fa t k and G^.m- 
This alternative is based on iteratively updating the matrix exponential from the computation 
at the previous step. We will first demonstrate its full version, and then provide a simple 
approximation which is not exactly symplectic on all variables but symplectic on the fast variables 
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(in the sense of a symplectic submanifold) and exhibits satisfactory long time performance in 
numerical experiments. 



Lemma 5.1. Define 



a(t) p(t) 7 (t) ' 
F 2 (t) G 2 (t) 
F 3 (t) 



exp 



-N T MJ ' 
-N T M 
ON 



(71) 



Then for any H, we have -F 3 (H) T j{H) = F 3 T (s)M(-JG 2 {t)) ds. 

Proof. Differentiating (|7ip with respect to t and equating each matrix component on left and 
right hand sides, we obtain: 

'a = -N T a 
F 2 = -N T F 2 
F 3 = NF 3 

f3 = -N T P + MJF 2 
G 2 = -N T G 2 + MF 3 
j = -N T 1 + MJG 2 



(72) 



where the initial conditions obviously are a(0) = I,F 2 (0) = I,F 3 (0) = I,f3(0) — 0,6*2(0) = 
0,7(0) = 0. 

Solving these inhomogeneous linear equations leads to known results including F 2 {t) — 
exp(-N T t) : F 3 (t) = exp(Nt) and G 2 {t) = J* exp(-N T (t - s))M exp(Ns) ds, as well as new 
results such as 



which is equivalent to 



7(t)= / expi-N 1 (t- s))MJG 2 {s)ds, 
Jo 



F 3 {H) T 1 {H)= F 3 (s) T M(-JG 2 (s))ds 
Jo 



(73) 



(74) 



□ 



Lemma 5.2. If M = M T , F^F 3 — I and dF 3 = —JG 2 , such as those derived from N and M 



defined in Section \2.3l then 



dG 2 (H) = F 2 (H) ^-{F 3 {H) T 1 (H)) T - F 3 (H) T -f(H) + ^ F 3 (sfdMF 3 {s) ds - {-JG 2 (H)) T G 2 (H) 

(75) 



Proof. By Leibniz's rule 

8G 2 {H) = [FsiHf]- 1 (d {F 3 {H) T G 2 (H))~dF 3 (H) T G 2 {H)) 
By the definition of F 3 and G 2l this is 

8G 2 (H) = F 2 (H) (d (jf F 3 {sfMF 3 (s) ds^j - {- J G 2 (H)) T G 2 {H)^j 



(76) 



(77) 
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in which 



/ pH \ pH pH pH 

d / F 3 {s) T MF 3 {s)ds\ = dF 3 (s) T MF 3 (s)ds + F 3 {sf MdF 3 (s) ds + F 3 (s) T dMF 3 {s) ds 
\Jo /Jo Jo Jo 



H 



H 



H 



for 7(i?) defined in Lemma |5. II 



{-JG 2 (s)) T MF 3 (s)ds + / F 3 (s) T M(-JG 2 {s))ds + / F 3 (s) T dMF 3 {s) ds 



(78) 
□ 



= -(F 3 (H) T 7 (H)) T - F 3 (if) T 7(ff) + / F 3 {sfdMF 3 (s) ds 

Jo 



Remark 5.1. J" F 3 (s) T dMF 3 {s) ds = F 3 (H) T G 2 (H) can be computed by again using the trick 
of: 



F 2 (t) G 2 (t) 
F 3 (t) 



exp 



-N T dM 
N 



(79) 



Of course, to get B' Q i j = djB§, l7 we use the fact that Bq^ = G 2 _o t i(H/2 n ). 

Lemma 5.3. Suppose q{™{ y , pg£ y , q?™ iy , pf™ x)l are obtained from q{? st \p{f st ,q^ w \pf}T 
by Integrator^ with Fs t k and G 2t k,i satisfying F^ k JF 3t k = J and G 2 _k.i — — J g g fL^ F 3 . k , then 



a slow v 
° g (fc+i)V' =I+ H 

pi ri slow o 



fast' 
fast 

LPfc' 



d 



G 2 ,k,j{H) JG 2 ,k,i{ H ) + F 3,k( H ) q siow G2 > k W 

oa k' 



fast' 
fast 



(80) 



Proof. Using chain rule, we have 
= I 



^slow 
i(k+l)',i 



H 
~2 



fast 
fast 



dqll 



This simplifies to (|8U)) because G 2: k 



— F 3 , k {HY G 2 ^{H) + F 3 ,k{HY 

k',j 



„slow 

3k', j 



' fast' 
fast 



-j 



F 3 , fc and —J 1 



J. 



(81) 
□ 



Integrator 3. Iterative matrix exponentiation scheme (alternative to Integrator^ that obtains 
F 3jk and G 2 ^,i via symplectic updates, k is the same index as the one used in Integrator]]^ n > 1 
is an integer controlling the accuracy of matrix exponential approximations. 



1. At the beginning of simulation, let q^, ow — q ilow - 
diK := Q^K(q s l ° w ) (i=l,..., d s ). Calculate 



-Hp s low and evaluate K a := K(q° l ° w ) and 



Co 



exp 



-No 1 M 0)S 
N a 



H/T 



by any favorite matrix exponentiation method (e.g., by the symplectic method introduced in 



Section [2~3)) , where Nq 



o I 

to 



and Mq ., 



e^diKo 



2. Compute B' i ■ = „ f low B$ j. One cheap way to do so is to use Lemma \5.2\ with Remark 



3. Start the updating loop, with the step count indicated by k starting from 1: let q^ low 'f ast 



% 



slow , fast slow, fast slow , fast 



Pi 



Po 



, and 



Sn — 



24 



4- Carry out the qk,Pk >-+ Qk',Pk' half-step (in Integrator^. Evaluate Kk '■= K(q k , ow ), and 



let D k 



e-^Kl-K^H/T 




Define Ak '■= A/.—1 exp(-Dfc) and use the equality 



cxp(Dk) = I + Dk (since Dk is nilpotent); similarly, define Ck '■= Ck-i exp(— D k ) 
Cfc-i — Ck-iDT ; 



5. Let Bk,i = — J g ®^iou, , which can be computed from known values using chain rule: 



B, 



-J 



a( c fc _i(j + p h )) 

f) n slow 

U( lk',i 

-slow 



AT. 



j=l U( ik',i U 1k-1' 

dD k 



{I + D k ) + C k -i 



dD k 
dq 



k' A 



• =1 u 1k'.i U( lk'. 



(82) 



To compute q^Iow > we need the derivatives of Kk and Kk—i with respect to q k ,° 
former is trivial, and the latter again can be computed by chain rule: 



the 



OKU _ * BtfSyj QKU 



~SlOW £)„SlOW 



(83) 



k' A 



3=1 



ifs',i 



ik-l',j 



6. B' k i ■ can be similarly computed from B' k _ 1 i ■ , Bk-\,i, C'k-i and Dk by repetitively applying 
chain rule. The detail is lengthy and hence omitted. 



7. Let k '■= Ak, G\ ki := Bk.i, F^ k := Ck, then repetitively apply 



713+1 rij+i 
2,) 




r 2i Lz 2,fc,i 



F 3+l 
r 3,k 



T?3 ri3 
r 2,k u 2,i,j 

H, k 



F 3 p3 p3 Q3 , Q3 p3 
r 2,k r 2.k r 2.k { - r 2,k,i ' { ~ r 2,k,i r 3,k 

F 3 F j 



for j 



. , n, and finally define 



^2,fe — r 2 k , Cr 2 ,fc,i .- Gr 2 M , ^ 3 , fe - t z k 



Compute 



by using Lemma \5.3l so that it could be used by Step\5\for the next k. 



—G2,k,i{H) is computed based on the following: 



-{AkBk.i + Bk^Ck) — —AkBk jJAkBk.i + AkB' ki 7 - + B' k ijCk — Bk^JBkj (84) 



where the first term is due to 



dA k 



= -AuBT . JA k , which is because dA T C + A T dC = 



d(A T C) = dl = and therefore dA T = ^A T dCC^ 1 = A T JBC~ X = A T JBA T . A similar 
trick of self multiplication applies to get the derivative of the 2" -times product. 

9. Carry out the qk',Pk' >-> qk+i,Pk+i half-step update of numerical integration using 1*2, fe ; 
Fs t k and Gi.k.i, and then increase k by 1 and go to Step^until integration time is reached. 



F^ ft and C?2,fe,i computed in this way (Integrator [3]) will also satisfy Lemma 13. II and render 
the integration symplectic on all variables. Proofs are omitted but they are analogous to those 
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in Section [XTJ and all structures, such as reversibility, symplecticity of F2 and F3 (illustrated by 
corresponding lemmas), and the relation between F3 and G2, will be preserved as long as they 
are satisfied by Aq, Bq^,C$ (i.e., the initial matrix exponentiation is accurate). 

In terms of efficiency, this method only uses one single matrix exponentiation operation and 
then keeps on updating it. Nevertheless, it is not easy to implement, and its speed advantage 
is not dominant. However, if the requirement on symplecticity is not that strict and a small 
numerical error in the matrix exponential is allowed (recall an analogous case of the famous 
implicit mid-point integrator, in which implicit solves are in fact not done perfectly and contin- 
uously polluting the symplecticity), we could use the approximation of q l k tilw'* = I- This will 

introduce a local error of 0(Hn/2 n ) in G2,k,i at each timestep (details omitted), but the local 
error in i<2,fc and F^ tk is 0, and the method is symplectic on the submanifold of the fast variables 
(although not symplectic on all variables). The approximating method is: 

Integrator 4. An efficient approximation of Integrator^ 

1. At the beginning of simulation, let q§ ow = q^ low rr^low 



diK := diK(q^ l , ow ) (i = 1, . . . , d s ) and calculate 



HpQ w and evaluate Kq 
A B 0< . 

exp 



Co 







by any favorite matrix exponentiation method, where Nq 



(I / 

~ l K 



K(q$ ow ) and 
and Mq i — 



e- l diK 6 



let q x 



slow , fast 



slow, fast 1 slow, fast 

q ana p ± 



slow, fast 

Po 





2. Start the updating loop, with the step count indicated by k starting from 1; 



3. Carry out the q k ,p k ► Ik' ,Pk' half-step. Evaluate K k 
"0 e-\Kl-Kl_ x )Hl^ 




K{q s k l , ow ) andd,K k := diK{qf}, 



slow \ 



let D k := 



and Ek 



e-^idiKk-diKk^H/T 1 



fine 

I + D k 


A k -iD k , 



~A k 


Bk,i 




Ak-i B k ^i_i 




D k 


E ky i 





Ck _ 







xexp 








(because D k E k j 



B k s — -Bfc-i,i + A k -\E k .i — B k -i yi DT , and C k 





and use the equality exp 

and E kt iD\ = 0) to evaluate A k 
Cfc-i — C k ~\DT; 




D k 




De- 

E k ,i 



A k -i 



Let Fl k 



A k, G 2ki :- B k ^, F 3;fc 



1 2 



C k , then repetitively apply 



r 2,k 



G 

F: 



i+i 

2.k.\ 



3,k 



F j 
r 2,k 






pi pi 

r 2,k r 2,k 


t?3 r<i _i_ r<3 wi 

r 2,k (Jr 2,k,i ' ^2,k,i r 3,k 





r 3,k _ 







pi pi 

r 3,k r 3,k 


F2,k '■= 


pn+l' 
r 2,k ' 


<-*2,fc,i ■— Lr 2 ,fe,i' 


^3,fc - -f 3ifc • 



for j = 1, . . . , n, and finally define 



5. Carry out the qk',Pk' <Ik+i,Pk+i half-step update of numerical integration using F2,k, 
F3 t k and G2,k,i, and then increase k by 1 and go to Step\^until integration time is reached. 

Numerical experiments presented in Section [4] are repeated using this approximating inte- 
grator. Energy preservations are as good as before, and slow trajectories show no significant 
deviation, suggesting no significant effect of the approximated symplecticity (detailed results 
omitted). This approximation, on the other hand, allows a choice of an even smaller n, such as 
n = 5 for the previous examples, which results in a further speed-up. 
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